{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "''' \n",
    "这里演示的是 不同预报方程对轨迹的影响。\n",
    "'''\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import solve_ivp\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "# 定义 Lorenz 63 系统的微分方程\n",
    "def lorenz63(t, state, sigma, rho, beta):\n",
    "    x, y, z = state\n",
    "    dx_dt = sigma * (y - x)\n",
    "    dy_dt = x * (rho - z) - y\n",
    "    dz_dt = x * y - beta * z\n",
    "    return [dx_dt, dy_dt, dz_dt]\n",
    "\n",
    "def lorenz63_fake(t, state, sigma, rho, beta):\n",
    "    \"\"\"假设参数估计错误\"\"\"\n",
    "    x, y, z = state\n",
    "    \n",
    "    # 假设我们对参数的估计有 5% 的误差\n",
    "    sigma_est = sigma * 1.05  # Prandtl 数估计偏高\n",
    "    rho_est = rho * 0.95      # Rayleigh 数估计偏低\n",
    "    beta_est = beta * 1.02    # 几何参数略有偏差\n",
    "    \n",
    "    dx_dt = sigma_est * (y - x)\n",
    "    dy_dt = x * (rho_est - z) - y\n",
    "    dz_dt = x * y - beta_est * z\n",
    "    \n",
    "    return [dx_dt, dy_dt, dz_dt]\n",
    "\n",
    "# 参数设置 (经典混沌参数)\n",
    "sigma = 10.0\n",
    "rho = 28.0\n",
    "beta = 8.0 / 3.0\n",
    "\n",
    "# 两个不同的初始条件\n",
    "initial_state1 = [1.0, 1.0, 1.0]                # 初始条件1\n",
    "initial_state2 = [1.0, 1.0, 1.0]        # 初始条件2 (微小差异)\n",
    "\n",
    "# 时间范围\n",
    "t_span = (0, 1)  # 缩短时间以更好观察差异\n",
    "t_eval = np.linspace(t_span[0], t_span[1], 3000)\n",
    "\n",
    "# 求解两个不同初始条件的微分方程\n",
    "solution1 = solve_ivp(lorenz63, t_span, initial_state1, args=(sigma, rho, beta), \n",
    "                      t_eval=t_eval, method='RK45', rtol=1e-6, atol=1e-9)\n",
    "\n",
    "solution2 = solve_ivp(lorenz63_fake, t_span, initial_state2, args=(sigma, rho, beta), \n",
    "                      t_eval=t_eval, method='RK45', rtol=1e-6, atol=1e-9)\n",
    "\n",
    "# 提取解\n",
    "x1, y1, z1 = solution1.y\n",
    "x2, y2, z2 = solution2.y\n",
    "\n",
    "# 创建图形\n",
    "plt.figure(figsize=(18, 7))\n",
    "\n",
    "# 图1: 3D 轨迹对比 - 合并在一个3D图中\n",
    "ax1 = plt.subplot(1, 2, 1, projection='3d')\n",
    "\n",
    "# 绘制两条轨迹，使用不同颜色\n",
    "ax1.plot(x1, y1, z1, lw=1.0, color='blue', alpha=0.8, label='Lorenz63')\n",
    "ax1.plot(x2, y2, z2, lw=1.0, color='red', alpha=0.8, label='lorenz63_fake')\n",
    "\n",
    "# 标记起始点\n",
    "ax1.scatter([initial_state1[0]], [initial_state1[1]], [initial_state1[2]], \n",
    "           color='blue', s=80, marker='o', label='Start Point 1', edgecolor='darkblue', linewidth=2)\n",
    "ax1.scatter([initial_state2[0]], [initial_state2[1]], [initial_state2[2]], \n",
    "           color='red', s=80, marker='s', label='Start Point 2', edgecolor='darkred', linewidth=2)\n",
    "\n",
    "# 标记结束点\n",
    "ax1.scatter([x1[-1]], [y1[-1]], [z1[-1]], \n",
    "           color='blue', s=100, marker='*', alpha=0.9, edgecolor='darkblue', linewidth=2)\n",
    "ax1.scatter([x2[-1]], [y2[-1]], [z2[-1]], \n",
    "           color='red', s=100, marker='*', alpha=0.9, edgecolor='darkred', linewidth=2)\n",
    "\n",
    "ax1.set_xlabel('X', fontsize=12)\n",
    "ax1.set_ylabel('Y', fontsize=12)\n",
    "ax1.set_zlabel('Z', fontsize=12)\n",
    "ax1.set_title('Lorenz 63 Attractor: Two Trajectories\\nInitial Difference: {:.2e}'.format(\n",
    "    abs(initial_state2[0] - initial_state1[0])), fontsize=14)\n",
    "ax1.legend(loc='upper right')\n",
    "\n",
    "# 设置视角以更好地观察差异\n",
    "ax1.view_init(elev=10, azim=45)\n",
    "\n",
    "# 图2: 蝴蝶效应 - 轨迹差异随时间的演化\n",
    "plt.subplot(1, 2, 2)\n",
    "\n",
    "# 计算两轨迹间的欧几里得距离\n",
    "distance = np.sqrt((x1 - x2)**2 + (y1 - y2)**2 + (z1 - z2)**2)\n",
    "\n",
    "# 绘制距离演化（对数尺度）\n",
    "plt.semilogy(solution1.t, distance, 'green', linewidth=2.5, label='Euclidean Distance')\n",
    "plt.xlabel('Time', fontsize=12)\n",
    "plt.ylabel('Distance (log scale)', fontsize=12)\n",
    "plt.title('Butterfly Effect: Trajectory Divergence', fontsize=14)\n",
    "plt.grid(True, alpha=0.3)\n",
    "plt.legend()\n",
    "\n",
    "# 添加参数信息\n",
    "textstr = f'Parameters:\\nσ = {sigma}\\nρ = {rho}\\nβ = {beta:.3f}\\n\\nInitial Conditions:\\nIC1: ({initial_state1[0]:.1f}, {initial_state1[1]:.1f}, {initial_state1[2]:.1f})\\nIC2: ({initial_state2[0]:.6f}, {initial_state2[1]:.1f}, {initial_state2[2]:.1f})'\n",
    "plt.text(0.02, 0.98, textstr, transform=plt.gca().transAxes, fontsize=10,\n",
    "         verticalalignment='top', bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.8))\n",
    "\n",
    "# # 添加一些关键点标注\n",
    "# max_distance_idx = np.argmax(distance)\n",
    "# plt.annotate(f'Max distance: {distance[max_distance_idx]:.3f}\\nat t = {solution1.t[max_distance_idx]:.1f}',\n",
    "#              xy=(solution1.t[max_distance_idx], distance[max_distance_idx]),\n",
    "#              xytext=(solution1.t[max_distance_idx] + 5, distance[max_distance_idx] * 0.1),\n",
    "#              arrowprops=dict(arrowstyle='->', color='red', alpha=0.7),\n",
    "#              bbox=dict(boxstyle='round,pad=0.3', facecolor='yellow', alpha=0.7))\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import solve_ivp\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "# 定义 Lorenz 63 系统的微分方程\n",
    "def lorenz63(t, state, sigma, rho, beta):\n",
    "    x, y, z = state\n",
    "    dx_dt = sigma * (y - x)\n",
    "    dy_dt = x * (rho - z) - y\n",
    "    dz_dt = x * y - beta * z\n",
    "    return [dx_dt, dy_dt, dz_dt]\n",
    "\n",
    "# 参数设置 (经典混沌参数)\n",
    "sigma = 10.0\n",
    "rho = 28.0\n",
    "beta = 8.0 / 3.0\n",
    "\n",
    "# 初始条件\n",
    "initial_state = [1.0, 1.0, 1.0]\n",
    "\n",
    "# 时间范围\n",
    "t_span = (0, 1)\n",
    "t_eval = np.linspace(t_span[0], t_span[1], 5000)\n",
    "\n",
    "# 求解微分方程\n",
    "solution = solve_ivp(lorenz63, t_span, initial_state, args=(sigma, rho, beta), \n",
    "                     t_eval=t_eval, method='RK45', rtol=1e-6, atol=1e-9)\n",
    "\n",
    "# 提取解\n",
    "x = solution.y[0]\n",
    "y = solution.y[1]\n",
    "z = solution.y[2]\n",
    "\n",
    "# 创建图形\n",
    "plt.figure(figsize=(15, 8))\n",
    "\n",
    "# 3D 轨迹图\n",
    "ax1 = plt.subplot(2, 2, 1, projection='3d')\n",
    "ax1.plot(x, y, z, lw=0.5, color='blue')\n",
    "ax1.set_xlabel('X')\n",
    "ax1.set_ylabel('Y')\n",
    "ax1.set_zlabel('Z')\n",
    "ax1.set_title('Lorenz 63 Attractor 3D View')\n",
    "\n",
    "# X-Y 平面投影\n",
    "plt.subplot(2, 2, 2)\n",
    "plt.plot(x, y, lw=0.5, color='red')\n",
    "plt.xlabel('X')\n",
    "plt.ylabel('Y')\n",
    "plt.title('X-Y Plane Projection')\n",
    "\n",
    "# 时间序列\n",
    "plt.subplot(2, 2, 3)\n",
    "plt.plot(solution.t, x, label='X', lw=0.5)\n",
    "plt.plot(solution.t, y, label='Y', lw=0.5)\n",
    "plt.plot(solution.t, z, label='Z', lw=0.5)\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('Value')\n",
    "plt.legend()\n",
    "plt.title('Time Series')\n",
    "\n",
    "# 对初始条件敏感性的演示 (蝴蝶效应)\n",
    "initial_state2 = [1.0 + 1e-5, 1.0, 1.0]  # 微小变化的初始条件\n",
    "solution2 = solve_ivp(lorenz63, t_span, initial_state2, args=(sigma, rho, beta), \n",
    "                      t_eval=t_eval, method='RK45', rtol=1e-6, atol=1e-9)\n",
    "\n",
    "plt.subplot(2, 2, 4)\n",
    "plt.plot(solution.t, x, label='Initial Condition 1', lw=1)\n",
    "plt.plot(solution2.t, solution2.y[0], label='Initial Condition 2', lw=1)\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('X Value')\n",
    "plt.legend()\n",
    "plt.title('Butterfly Effect: Sensitivity to Initial Conditions')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import solve_ivp\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "# 定义 Lorenz 63 系统的微分方程\n",
    "def lorenz63(t, state, sigma, rho, beta):\n",
    "    x, y, z = state\n",
    "    dx_dt = sigma * (y - x)\n",
    "    dy_dt = x * (rho - z) - y\n",
    "    dz_dt = x * y - beta * z\n",
    "    return [dx_dt, dy_dt, dz_dt]\n",
    "\n",
    "# 参数设置 (经典混沌参数)\n",
    "sigma = 10.0\n",
    "rho = 28.0\n",
    "beta = 8.0 / 3.0\n",
    "\n",
    "# 初始条件\n",
    "initial_state1 = [1.0, 1.0, 1.0]\n",
    "initial_state2 = [1.0 + 1e-5, 1.0, 1.0]  # 微小变化的初始条件\n",
    "\n",
    "# 时间范围\n",
    "t_span = (0, 50)\n",
    "t_eval = np.linspace(t_span[0], t_span[1], 5000)\n",
    "\n",
    "# 求解微分方程 - 第一个初始条件\n",
    "solution1 = solve_ivp(lorenz63, t_span, initial_state1, args=(sigma, rho, beta), \n",
    "                     t_eval=t_eval, method='RK45', rtol=1e-6, atol=1e-9)\n",
    "\n",
    "# 求解微分方程 - 第二个初始条件\n",
    "solution2 = solve_ivp(lorenz63, t_span, initial_state2, args=(sigma, rho, beta), \n",
    "                     t_eval=t_eval, method='RK45', rtol=1e-6, atol=1e-9)\n",
    "\n",
    "# 提取解\n",
    "x1, y1, z1 = solution1.y\n",
    "x2, y2, z2 = solution2.y\n",
    "\n",
    "# 创建图形\n",
    "plt.figure(figsize=(15, 12))\n",
    "\n",
    "# 第一行：两个3D轨迹图\n",
    "# 第一个初始条件的3D轨迹\n",
    "ax1 = plt.subplot(3, 2, 1, projection='3d')\n",
    "ax1.plot(x1, y1, z1, lw=0.5, color='blue')\n",
    "ax1.set_xlabel('X')\n",
    "ax1.set_ylabel('Y')\n",
    "ax1.set_zlabel('Z')\n",
    "ax1.set_title('Initial Condition 1: [1.0, 1.0, 1.0]')\n",
    "\n",
    "# 第二个初始条件的3D轨迹\n",
    "ax2 = plt.subplot(3, 2, 2, projection='3d')\n",
    "ax2.plot(x2, y2, z2, lw=0.5, color='red')\n",
    "ax2.set_xlabel('X')\n",
    "ax2.set_ylabel('Y')\n",
    "ax2.set_zlabel('Z')\n",
    "ax2.set_title('Initial Condition 2: [1.00001, 1.0, 1.0]')\n",
    "\n",
    "# 第二行：X维度的时间序列比较\n",
    "plt.subplot(3, 2, 3)\n",
    "plt.plot(solution1.t, x1, label='IC1: X', lw=1, color='blue')\n",
    "plt.plot(solution2.t, x2, label='IC2: X', lw=1, color='red')\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('X Value')\n",
    "plt.legend()\n",
    "plt.title('X Dimension Comparison')\n",
    "plt.grid(True, alpha=0.3)\n",
    "\n",
    "# 第三行：Y维度的时间序列比较\n",
    "plt.subplot(3, 2, 4)\n",
    "plt.plot(solution1.t, y1, label='IC1: Y', lw=1, color='blue')\n",
    "plt.plot(solution2.t, y2, label='IC2: Y', lw=1, color='red')\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('Y Value')\n",
    "plt.legend()\n",
    "plt.title('Y Dimension Comparison')\n",
    "plt.grid(True, alpha=0.3)\n",
    "\n",
    "# 第四行：Z维度的时间序列比较\n",
    "plt.subplot(3, 2, 5)\n",
    "plt.plot(solution1.t, z1, label='IC1: Z', lw=1, color='blue')\n",
    "plt.plot(solution2.t, z2, label='IC2: Z', lw=1, color='red')\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('Z Value')\n",
    "plt.legend()\n",
    "plt.title('Z Dimension Comparison')\n",
    "plt.grid(True, alpha=0.3)\n",
    "\n",
    "# 第六个子图：显示两个轨迹的差异（绝对值）\n",
    "plt.subplot(3, 2, 6)\n",
    "diff_x = np.abs(x1 - x2)\n",
    "diff_y = np.abs(y1 - y2)\n",
    "diff_z = np.abs(z1 - z2)\n",
    "plt.plot(solution1.t, diff_x, label='|X1 - X2|', lw=1, alpha=0.8)\n",
    "plt.plot(solution1.t, diff_y, label='|Y1 - Y2|', lw=1, alpha=0.8)\n",
    "plt.plot(solution1.t, diff_z, label='|Z1 - Z2|', lw=1, alpha=0.8)\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('Absolute Difference')\n",
    "plt.legend()\n",
    "plt.title('Absolute Differences Between Trajectories')\n",
    "plt.yscale('log')  # 使用对数坐标更好地显示指数增长\n",
    "plt.grid(True, alpha=0.3)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "gsj",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
